Exercise 2: Can you improve these priors?
Exercise 4: Coding the logistic regression
Exercise 5: Coding the logistic regression to run in parallel (optional)
Exercise 6: Summarizing coda objects
Exercise 7: Understanding coda objects
Exercise 8: Convert the zm object to a data frame
Exercise 9: Vectors in coda objects
Exercise 10: Making plots with JAGS objects
Exercise 11: Summarizing the JAGS object
Exercise 12: Assessing convergence
You can complete all the coding exercises in the JAGS Primer using simple R scripts. However, you might be wondering about how we created the JAGS Primer key you are looking at right now. We did this using Yihui Xie’s knitr package in R. This can be a highly useful tools for organizing your Bayesian analyses. Within the same document, you can:
Best of all, knitr can produce beautiful html files as output, which can be easily shared with collaborators. We encourage you to become familiar with knitr. We reccomend Karl Broman’s knitr in a nutshell as an excellent introductory tutorial. You can also open JagsPrimerAnswers.Rmd to see how this html document is generated.
Why does \(x\) fail to appear in the posterior distribution? Draw the Bayesian network for this model.
Answer: There is no \(x\) because we are assuming it is measured without error.
A recurring theme in this course will be to use priors that are informative whenever possible. The gamma priors in equation 3 include the entire number line \(>0\). Don’t we know more about population biology than that? Lets, say for now that we are modeling the population dynamics of a large mammal. How might you go about making the priors on population parameters more informative?
Answer: A great source for priors in biology are allometric scaling relationships that predict all kinds of biological quantities based on the mass organisms (Peters, 1983; Pennycuick,1992). If you know the approximate mass of the animal, you can compute broad but nonetheless informative priors on \(r\) and \(K\). This might leave out the social scientists, but I would trust the scaling of \(r\) for people if not \(K\).
In the absence of some sort of scholarly way to find priors, we can at least constrain them somewhat. There is no way that a large mammal can have an intrinsic rate of increase exceeding 1 – many values for \(r\) within gamma(.001, .001) are far large than than that and hence are complete nonsense. We know \(r\) must be positive and we can put a plausible bound on its upper limit. The only requirement for a vague prior is that its “\(\ldots\) range of uncertainty should be clearly wider that the range of reasonable values of the parameter\(\ldots\)” (Gelman and Hill, 2009, page 355), so we could use \(r\) ~ uniform(0, 2) and be sure that it would be minimally informative. Similarly, we could use experience and knowledge to put some reasonable bounds on \(K\) and even \(\sigma\), which we can use to calculate \(\tau\) as \(\tau=\frac{1}{\sigma^{2}}\).
Peters. The ecological implications of body size. Cambridge University Press, Cambridge, United Kingdom, 1983.
C. J. Pennycuick. Newton rules biology. Oxford University Press, Oxford United Kingdom, 1992.
A. Gelman and J. Hill. Data analysis using regression and multilievel / hierarchical modeling. Cambridge University Press, Cambridge, United Kingdom, 2009.
for loopsWrite a code fragment to set vague normal priors for 5 regression coefficients – dnorm(0, 10E-6) – stored in the vector b.
for(i in 1:5){
b[i] ~ dnorm(0, .000001)
}
Write R code (algorithm 3) to run the JAGS model (algorithm 2) and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). We suggest you insert the JAGS model into this R script using the sink command as shown in algorithm 4. You will find this a very convenient way to keep all your code in the same R script.
Here is the joint distribution for our logisitic model again, with the priors updated from exercise 2 and \(\tau\) expressed as a derived quantity,
\[\begin{eqnarray} \mu_{i} & = & r-\frac{rx_{i}}{K}\textrm{,}\nonumber\\[1em] \tau & = & \frac{1}{\sigma^{2}}\textrm{,}\nonumber\\[1em] \left[r,K,\sigma\mid\mathbf{y}\right] & \propto & \prod_{i=1}^{n}\textrm{normal}\left(y_{i} \mid \mu_{i},\tau\right)\textrm{uniform}\left(K\mid 0,4000\right) \textrm{uniform}\left(\sigma\mid 0, 2\right) \textrm{uniform}\left(r\mid 0,2\right)\textrm{.}\nonumber\\ \end{eqnarray}\]We use the sink command to create a JAGS script from our joint distribution. This file is created within R and saved in the working directory. Please note that the outer set of brackets are only required when running this code within an R markdown document (as we did to make this answer key). If you are running them in a plain R script, they are not needed.
It might be useful for you to know that the data were generated from a model with \(r=.2\), \(K=1200\) and \(\sigma = .03\)
{ # Extra bracket needed only for R markdown files
sink("LogisticJAGS.R")
cat("
model{
# priors
K ~ dunif(0, 4000)
r ~ dunif (0, 2)
sigma ~ dunif(0, 2)
tau <- 1/sigma^2
# likelihood
for(i in 1:n){
mu[i] <- r - r/K * x[i]
y[i] ~ dnorm(mu[i], tau)
}
}
",fill = TRUE)
sink()
} # Extra bracket needed only for R markdown files
Then we run the remaining commands discussed in the JAGS Primer. Note that jm is calling the JAGS script LogisticJAGS.R we just created.
rm(list = ls())
library(ESS575)
library(rjags)
## Loading required package: coda
## Linked to JAGS 4.2.0
## Loaded modules: basemod,bugs
set.seed(1)
Logistic <- Logistic[order(Logistic$PopulationSize),]
inits = list(
list(K = 1500, r = .2, sigma = 1),
list(K = 1000, r = .15, sigma = .1),
list(K = 900, r = .3, sigma = .01))
data = list(
n = nrow(Logistic),
x = as.double(Logistic$PopulationSize),
y = as.double(Logistic$GrowthRate))
n.adapt = 5000
n.update = 10000
n.iter = 10000
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 217
##
## Initializing model
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
summary(zm)
##
## Iterations = 15001:25000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## K 1.238e+03 6.288e+01 3.630e-01 7.539e-01
## r 2.006e-01 9.773e-03 5.643e-05 1.143e-04
## sigma 2.867e-02 3.033e-03 1.751e-05 2.522e-05
## tau 1.257e+03 2.594e+02 1.498e+00 2.042e+00
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## K 1.131e+03 1.195e+03 1233.0596 1.275e+03 1375.9685
## r 1.814e-01 1.941e-01 0.2006 2.071e-01 0.2199
## sigma 2.345e-02 2.654e-02 0.0284 3.052e-02 0.0354
## tau 7.980e+02 1.073e+03 1239.4604 1.420e+03 1819.0041
Append R code (algorithm 5) to the script you made in exercise 4 to run the JAGS model (algorithm 2) in parallel and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). Use the proc.time function in R to compare the time required for the sequential and parallel JAGS run. If your computer has 3 cores, try running only 2 chains in parallel when doing this exercise. If you have fewer than 3 cores, work with a classmate that has at least 3 cores.
We create a function called initFunc to randomly draw values from a portion of each parameter’s support. We will use this function to provide initial values for each chain. We test out initFunc by running it a couple of times.
initFunc <- function (){
return(list(
K = runif(1, 10, 4000),
r = runif(1, .01, 2),
sigma = runif(1, 0, 2)))}
initFunc()
## $K
## [1] 1069.38
##
## $r
## [1] 0.7505266
##
## $sigma
## [1] 1.145707
initFunc()
## $K
## [1] 3633.749
##
## $r
## [1] 0.411347
##
## $sigma
## [1] 1.796779
Now we run the model in parallel using the code from algorithm 5. We use the proc.time function to see how long JAGS takes to run the Logistics model in parallel.
# run JAGS model in parallel
library(parallel)
detectCores()
## [1] 8
cl <- makeCluster(3) # Here we use three cores
clusterExport(cl, c("data", "initFunc", "n.adapt", "n.update", "n.iter"))
ptm <- proc.time()
out <- clusterEvalQ(cl, {
library(rjags)
set.seed(1)
jm = jags.model("LogisticJAGS.R", data = data, inits = initFunc(),
n.chains = 1, n.adapt = n.adapt)
update(jm, n.iter = n.update)
zmCore = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zmCore))
})
ParallelTime <- proc.time() - ptm
ParallelTime
## user system elapsed
## 0.003 0.000 1.242
stopCluster(cl)
zmP <- mcmc.list(out)
We rerun the model sequentially and use proc.time again for comparison.
ptm <- proc.time()
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 217
##
## Initializing model
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
SequentialTime <- proc.time() - ptm
SequentialTime
## user system elapsed
## 4.144 0.008 4.154
Looks like the parallel model runs 3.34 times faster. This factor should increase the more iterations you run (why?) to a limit of 3.
Build a table that contains the mean, standard deviation, median and upper and lower 2.5% CI for the parameters from the logistic model. Output your table with 3 significant digits to .csv file readable by Excel (Hint: see the signif() function).
a <- signif(as.data.frame(summary(zm)$stat[, 1:2]), digits = 3)
b <- signif(as.data.frame(summary(zm)$quantile[, c(1, 3, 5)]), digits = 3)
LogisticParameters <- cbind(rownames(a), a, b)
rownames(LogisticParameters) <- c()
names(LogisticParameters) <- c("parameter", "mean", "standardDeviation", "lower95", "median", "upper95")
LogisticParameters
## parameter mean standardDeviation lower95 median upper95
## 1 K 1.24e+03 6.31e+01 1.13e+03 1.23e+03 1.38e+03
## 2 r 2.01e-01 9.76e-03 1.82e-01 2.01e-01 2.20e-01
## 3 sigma 2.87e-02 3.04e-03 2.35e-02 2.84e-02 3.54e-02
## 4 tau 1.25e+03 2.59e+02 7.99e+02 1.24e+03 1.81e+03
write.csv(LogisticParameters, file = "LogisticParameters.csv")
Modify your code to produce a coda object with 3 chains called zm.short, setting n.adapt = 500, n.update = 500, and n.iter = 20.
zm.short at the console.n.adapt = 500
n.update = 500
n.iter = 20
jm.short = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 217
##
## Initializing model
update(jm.short, n.iter = n.update)
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
zm.short[[2]][3,3] # third iteration from second chain for sigma
## sigma
## 0.03169468
zm.short[[1]][,2] # all estimates of r from first chain
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1001
## End = 1020
## Thinning interval = 1
## [1] 0.2030203 0.2108610 0.2202984 0.2189606 0.2078010 0.2008283 0.1961416
## [8] 0.1959804 0.2015194 0.2188144 0.2048849 0.1867363 0.1925875 0.1955791
## [15] 0.1957873 0.2034632 0.1941668 0.2004273 0.1993938 0.1917201
zm.short
## [[1]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1001
## End = 1020
## Thinning interval = 1
## K r sigma tau
## [1,] 1181.653 0.2030203 0.02404521 1729.5884
## [2,] 1235.919 0.2108610 0.03156689 1003.5439
## [3,] 1209.245 0.2202984 0.02936133 1159.9746
## [4,] 1173.462 0.2189606 0.02973226 1131.2126
## [5,] 1186.633 0.2078010 0.03252151 945.4937
## [6,] 1191.945 0.2008283 0.03339594 896.6286
## [7,] 1203.969 0.1961416 0.03030680 1088.7291
## [8,] 1191.896 0.1959804 0.02979303 1126.6025
## [9,] 1201.578 0.2015194 0.02506923 1591.1756
## [10,] 1197.818 0.2188144 0.02979680 1126.3176
## [11,] 1202.582 0.2048849 0.02979170 1126.7028
## [12,] 1290.592 0.1867363 0.02691262 1380.6636
## [13,] 1295.286 0.1925875 0.02779452 1294.4389
## [14,] 1285.023 0.1955791 0.02784065 1290.1529
## [15,] 1263.384 0.1957873 0.02681210 1391.0356
## [16,] 1304.803 0.2034632 0.02567445 1517.0426
## [17,] 1223.385 0.1941668 0.03049932 1075.0277
## [18,] 1238.984 0.2004273 0.03178119 990.0555
## [19,] 1210.264 0.1993938 0.02703035 1368.6629
## [20,] 1313.790 0.1917201 0.02513456 1582.9140
##
## [[2]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1001
## End = 1020
## Thinning interval = 1
## K r sigma tau
## [1,] 1259.236 0.1987620 0.02934469 1161.2906
## [2,] 1296.763 0.1782903 0.02826893 1251.3568
## [3,] 1361.928 0.1808955 0.03169468 995.4679
## [4,] 1341.617 0.1928334 0.02663881 1409.1928
## [5,] 1329.236 0.1800192 0.02653763 1419.9583
## [6,] 1307.419 0.1854929 0.02474384 1633.3000
## [7,] 1329.466 0.1883291 0.02707302 1364.3525
## [8,] 1260.794 0.2078104 0.02699141 1372.6155
## [9,] 1215.190 0.1978350 0.02504065 1594.8100
## [10,] 1252.989 0.2026278 0.02612994 1464.6136
## [11,] 1190.874 0.2154729 0.02303833 1884.0747
## [12,] 1164.754 0.1959241 0.02600548 1478.6662
## [13,] 1272.806 0.1976314 0.02665192 1407.8067
## [14,] 1197.249 0.2104054 0.02854932 1226.8977
## [15,] 1186.206 0.2095827 0.02748301 1323.9499
## [16,] 1197.686 0.2062743 0.03478085 826.6463
## [17,] 1232.439 0.2067848 0.03155640 1004.2112
## [18,] 1233.755 0.2092772 0.02854276 1227.4624
## [19,] 1260.815 0.2016824 0.02396908 1740.5932
## [20,] 1256.279 0.2028321 0.03251925 945.6253
##
## [[3]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1001
## End = 1020
## Thinning interval = 1
## K r sigma tau
## [1,] 1226.462 0.2014688 0.02563643 1521.5452
## [2,] 1223.745 0.2005959 0.02878547 1206.8500
## [3,] 1239.955 0.1964147 0.03268001 936.3442
## [4,] 1287.268 0.1849378 0.03336379 898.3573
## [5,] 1285.660 0.1837527 0.03266691 937.0957
## [6,] 1375.630 0.1927980 0.03318887 907.8522
## [7,] 1383.053 0.1927733 0.03195796 979.1333
## [8,] 1356.639 0.1934964 0.02987360 1120.5332
## [9,] 1319.808 0.1904076 0.03017143 1098.5208
## [10,] 1281.788 0.1885839 0.02713142 1358.4854
## [11,] 1398.285 0.1864175 0.02899040 1189.8485
## [12,] 1320.734 0.1965401 0.02913553 1178.0244
## [13,] 1231.152 0.1949804 0.03286307 925.9420
## [14,] 1256.216 0.1892845 0.03561903 788.1987
## [15,] 1332.042 0.1909689 0.02846402 1234.2621
## [16,] 1322.767 0.1905138 0.02875841 1209.1224
## [17,] 1279.100 0.1871464 0.02947830 1150.7872
## [18,] 1370.564 0.1964373 0.03139866 1014.3267
## [19,] 1245.975 0.1987005 0.03006883 1106.0303
## [20,] 1165.587 0.1989103 0.02571081 1512.7543
##
## attr(,"class")
## [1] "mcmc.list"
zm object to a data frameUsing the elements of data frame (not zm) as input to functions:
ecdf function.)df = as.data.frame(rbind(zm[[1]], zm[[2]], zm[[3]]))
# Find the maximum value of sigma
max(df$sigma)
## [1] 0.04593981
# Find the mean of r for the first 1000 iterations
mean(df$r[1: 1000])
## [1] 0.20099
# Find the mean of r for the first 1000 iterations
nr = length(df$r)
mean(df$r[(nr - 1000): nr])
## [1] 0.2011763
plot(density(df$K), main = "", xlim=c(800, 2000), xlab = "K")
Fig. 2. Posterior density of K.
# Find the probability that the parameter K exceeds 1600
1 - ecdf(df$K)(1600)
## [1] 0.0001333333
# Find the probability that the parameter 1000 < K < 1300
ecdf(df$K)(1300) - ecdf(df$K)(1000)
## [1] 0.8486667
Modify your code to include estimates of \(\mu\) and summarize the coda object zm. What if you wanted to plot the model predictions with 95% credible intervals against the data. How would you do that?
Answer: There are several ways this can be done, but the general idea is that you need to extract the rows of the coda object that contain the quantiles for \(\mu\), which can be tedious and error prone. For example, if you use rows in the summary table and add or subtract parameters to be estimated, then your row counts will be off. There are ways to use rownames, but a far better way to plot vectors is described in the section on JAGS objects.
n.adapt = 5000
n.update = 10000
n.iter = 10000
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 217
##
## Initializing model
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "mu"), n.iter = n.iter, n.thin = 1)
zj = jags.samples(jm, variable.names = c("K", "r", "sigma", "mu"), n.iter = n.iter, n.thin = 1)
mu <- as.data.frame(summary(zm)$quantile[2:51, c(1, 3, 5)]) # <- this is an easy place to make a mistake
names(mu) <- c("lower95", "median", "upper95")
Logistic2 <- cbind(mu, Logistic)
plot(Logistic2$PopulationSize, Logistic2$GrowthRate, xlab = "N", ylab = "Per capita growth rate")
lines(Logistic2$PopulationSize, Logistic2$median)
lines(Logistic2$PopulationSize, Logistic2$lower95, lty = "dashed")
lines(Logistic2$PopulationSize, Logistic2$upper95, lty = "dashed")
Fig. 3. Median and 95% credible intervals for predicted growth rate.
For the logistic model:
For convenience, we created the JAGS object zj in exercise 9. Here we use it to make plots.
mu <- summary(zj$mu, quantile, c(.025, .5, .975))$stat # <- notice this is harder to mess up!
par(mfrow = c(1, 2))
plot(Logistic$PopulationSize, Logistic$GrowthRate, xlab = "N", ylab = "Per capita growth rate")
lines(Logistic$PopulationSize, mu[2,])
lines(Logistic$PopulationSize, mu[1,], lty = "dashed")
lines(Logistic$PopulationSize, mu[3,], lty = "dashed")
plot(density(zj$K), main = "", xlim=c(800, 2000), xlab = "K")
Fig. 4. Median and 95% credible intervals for predicted growth rate and posterior density of K.
summary function.summary(zj$K, median)$stat
## [1] 1233.173
quantile(zj$mu[16,,], c(.025, .975))
## 2.5% 97.5%
## 0.1179175 0.1362593
ecdf(zj$mu[16,,])(0)
## [1] 0
Rerun the logistic model with n.adapt = 100. Then do the following:
traceplot and with the Gelman-Rubin, Heidelberger and Welch, and Raftery diagnostics.set.seed(1)
n.adapt = 100
jm.short = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 217
##
## Initializing model
n.iter = 100
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
plot(zm.short)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.06 1.21
## r 1.10 1.30
## sigma 1.02 1.05
## tau 1.01 1.03
##
## Multivariate psrf
##
## 1.08
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.1290
## r failed NA 0.0401
## sigma passed 1 0.9767
## tau passed 1 0.9490
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.26e+03 3.03e+01
## r <NA> NA NA
## sigma passed 2.82e-02 7.42e-04
## tau passed 1.30e+03 5.73e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.439
## r passed 1 0.921
## sigma passed 1 0.451
## tau passed 1 0.513
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 5.02e+01
## r passed 2.01e-01 5.17e-03
## sigma passed 2.89e-02 7.96e-04
## tau passed 1.24e+03 7.05e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K failed NA 0.00132
## r failed NA 0.03980
## sigma passed 11 0.06082
## tau passed 1 0.05201
##
## Halfwidth Mean Halfwidth
## test
## K <NA> NA NA
## r <NA> NA NA
## sigma passed 2.81e-02 6.27e-04
## tau passed 1.28e+03 6.54e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
n.update = 500
update(jm.short, n.iter = n.update)
n.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
plot(zm.short)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.02 1.03
## r 1.01 1.01
## sigma 1.00 1.02
## tau 1.00 1.01
##
## Multivariate psrf
##
## 1.01
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.703
## r passed 1 0.541
## sigma passed 1 0.354
## tau passed 1 0.343
##
## Halfwidth Mean Halfwidth
## test
## K passed 1244.883 1.78e+01
## r passed 0.200 2.14e-03
## sigma passed 0.029 4.74e-04
## tau passed 1233.463 3.96e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 51 0.1310
## r passed 51 0.0892
## sigma passed 1 0.7220
## tau passed 1 0.8296
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 2.09e+01
## r passed 2.01e-01 2.70e-03
## sigma passed 2.86e-02 3.86e-04
## tau passed 1.26e+03 3.25e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.243
## r passed 1 0.599
## sigma passed 1 0.808
## tau passed 1 0.874
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.08e+01
## r passed 2.01e-01 1.61e-03
## sigma passed 2.86e-02 3.24e-04
## tau passed 1.26e+03 2.75e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
n.update = 500
update(jm.short, n.iter = n.update)
n.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
plot(zm.short)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.00 1.01
## r 1.01 1.03
## sigma 1.00 1.00
## tau 1.00 1.00
##
## Multivariate psrf
##
## 1.01
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.529
## r passed 1 0.518
## sigma passed 1 0.659
## tau passed 1 0.326
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 10.46395
## r passed 2.02e-01 0.00143
## sigma passed 2.86e-02 0.00034
## tau passed 1.27e+03 27.42649
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.998
## r passed 1 0.950
## sigma passed 1 0.559
## tau passed 1 0.553
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 1.59e+01
## r passed 2.01e-01 1.56e-03
## sigma passed 2.85e-02 3.24e-04
## tau passed 1.27e+03 2.74e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.294
## r passed 1 0.257
## sigma passed 51 0.140
## tau passed 1 0.139
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 1.15e+01
## r passed 2.00e-01 1.71e-03
## sigma passed 2.84e-02 3.35e-04
## tau passed 1.27e+03 2.68e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
#Final run--you know this has converged!
n.update = 5000
update(jm.short, n.iter = n.update)
n.iter = 10000
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
plot(zm.short)
gelman.diag(zm.short)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1 1
## r 1 1
## sigma 1 1
## tau 1 1
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.264
## r passed 1 0.153
## sigma passed 1 0.124
## tau passed 1 0.113
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 3.25e+00
## r passed 2.01e-01 4.35e-04
## sigma passed 2.86e-02 7.76e-05
## tau passed 1.26e+03 6.68e+00
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.641
## r passed 1 0.527
## sigma passed 1 0.294
## tau passed 1 0.186
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 4.81e+00
## r passed 2.01e-01 5.78e-04
## sigma passed 2.87e-02 8.05e-05
## tau passed 1.25e+03 6.59e+00
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.726
## r passed 1 0.833
## sigma passed 1 0.643
## tau passed 1 0.626
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 2.88e+00
## r passed 2.01e-01 3.72e-04
## sigma passed 2.87e-02 7.98e-05
## tau passed 1.26e+03 6.35e+00
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 11820 3746 3.16
## r 10 10984 3746 2.93
## sigma 4 4752 3746 1.27
## tau 5 5973 3746 1.59
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 10910 3746 2.91
## r 12 15558 3746 4.15
## sigma 4 5123 3746 1.37
## tau 6 6816 3746 1.82
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 8 8518 3746 2.27
## r 8 9508 3746 2.54
## sigma 4 5210 3746 1.39
## tau 6 7002 3746 1.87